/* 
 * File:   meshimpl_surf3d.cc
 * Author: boyan
 * 
 * Created on March 27, 2009, 4:26 PM
 */

#include "Mesh_grid3d.h"


#include "Mesh_surf3d.h"


#include <assert.h>
#include <cstdlib>
#include <iostream>
#include <algorithm>
#include <vector>

void Mesh_surf3d::surf_init(const std::vector<int> &in_elems)
{
    std::vector<int> foo;

    // make surface faces
    surf_faces.resize(my_grid->NF, false);
    foo.clear();
    foo.resize(my_grid->NF, 0);
    for (size_t idx = 0; idx < in_elems.size(); idx++)
    {
        int Tid = in_elems[idx];
        elem_t *T = my_grid->elems + Tid;
        for (int j = 0; j < 4; j++)
            ++foo[T->faces[j]];
    }
    for (size_t idx = 0; idx < foo.size(); idx++)
        switch (foo.at(idx)) // the number of elements containing face idx
        {
        case 1:
            surf_faces[idx] = true;
            break;
        case 0: // face is outside the patch of elements in in_elems
        case 2: // face is interior to the patch
            // in both case we leave the defauls `false`
            break;
        default:
            std::cerr << __FILE__ << ':' << __LINE__
                    << ": error making surface faces\n";
            exit(EXIT_FAILURE);
        }

    Mesquite::MsqError err;
    Mesquite::VertexIterator *f_iter = element_iterator(err); // must delete

    // decide normal signs
    swap.resize(surf_faces.size(), false);
    for (f_iter->restart(); !f_iter->is_at_end(); f_iter->operator ++())
    {
        int fid = handle2elem(f_iter->operator *());
        face_t *f = my_grid->faces + fid;
        std::vector<int>::const_iterator found;
        found = std::find(in_elems.begin(), in_elems.end(), f->elems[0]);
        if (found == in_elems.end())
            found = std::find(in_elems.begin(), in_elems.end(), f->elems[1]);
        if (found == in_elems.end())
        {
            std::cerr << __FILE__ << ':' << __LINE__
                    << ": error making normal signs\n";
            exit(EXIT_FAILURE);
        }
        int Tid = found.operator *();
        elem_t *T = my_grid->elems + Tid;
        for (int i = 0; i < 4; i++)
        {
            if (T->faces[i] == fid)
            {
                swap[fid] = elem_face_normal_sign(my_grid, T, i) < 0;
                break;
            }
        }
    }

    // make surface nodes
    surf_nodes.resize(my_grid->NN, false);
    foo.clear();
    foo.resize(my_grid->NN, 0);
    for (f_iter->restart(); !f_iter->is_at_end(); f_iter->operator ++())
    {
        int fid = handle2elem(f_iter->operator *());
        face_t *f = my_grid->faces + fid;
        for (size_t j = 0; j < 3; j++)
            ++foo[f->nodes[j]];
    }
    for (size_t idx = 0; idx < foo.size(); idx++)
        surf_nodes[idx] = foo.at(idx) > 0;

    delete f_iter;
}

Mesh_surf3d::Mesh_surf3d(grid3d_t &grid)
: Mesh_grid3d(grid)
{
    std::vector<int> in_el;
    in_el.reserve(my_grid->NEL);
    for (int i = 0; i < my_grid->NEL; i++)
        in_el.push_back(i);
    surf_init(in_el);
}

void Mesh_surf3d::get_all_elements(msq_std::vector<ElementHandle>& elements,
                                   Mesquite::MsqError& err)
{
    Mesquite::ElementIterator *it = element_iterator(err);
    int NFZ = 0;
    for (it->restart(); !it->is_at_end(); it->operator ++())
        ++NFZ;
    elements.clear();
    elements.reserve(NFZ);
    for (it->restart(); !it->is_at_end(); it->operator ++())
        elements.push_back(it->operator *());
    delete it;
}

void Mesh_surf3d::get_all_vertices(msq_std::vector<VertexHandle>& vertices,
                                   Mesquite::MsqError& err)
{
    Mesquite::VertexIterator *it = vertex_iterator(err);
    int NVZ = 0;
    for (it->restart(); !it->is_at_end(); it->operator ++())
        ++NVZ;
    vertices.clear();
    vertices.reserve(NVZ);
    for (it->restart(); !it->is_at_end(); it->operator ++())
        vertices.push_back(it->operator *());
    delete it;
}

void Mesh_surf3d::vertices_get_attached_elements
(
 const VertexHandle* vert_array, size_t num_vertex,
 msq_std::vector<ElementHandle>& elements,
 msq_std::vector<size_t>& offsets, Mesquite::MsqError& err)
{
    if (my_grid->inverse == NULL)
        MakeInverse(my_grid);
    elements.clear();
    offsets.clear();
    offsets.resize(num_vertex + 1);
    size_t current_offset = 0;
    for (size_t idx = 0; idx < num_vertex; idx++)
    {
        int surf_nid = handle2node(vert_array[idx]);
        assert(surf_nid < (int) surf_nodes.size());
        bool nid = surf_nodes.at(surf_nid);
        assert(nid == true);
        offsets[idx] = current_offset;
        node_inv_t *ninv = my_grid->inverse->ninv + surf_nid;
        for (int qq = 0; qq < ninv->NF; ++qq)
            if (is_face_on_surface(ninv->fids[qq]))
                ++current_offset;
    }
    offsets[num_vertex] = current_offset;
    elements.reserve(current_offset);
    for (size_t idx = 0; idx < num_vertex; idx++)
    {
        int surf_nid = handle2node(vert_array[idx]);
        node_inv_t *ninv = my_grid->inverse->ninv + surf_nid;
        for (int qq = 0; qq < ninv->NF; ++qq)
            if (is_face_on_surface(ninv->fids[qq]))
                elements.push_back(elem2handle(ninv->fids[qq]));
    }
    return;
}

void Mesh_surf3d::elements_get_attached_vertices
(
 const ElementHandle *elem_handles, size_t num_elems,
 msq_std::vector<VertexHandle>& vert_handles,
 msq_std::vector<size_t>& offsets, Mesquite::MsqError &err)
{
    vert_handles.clear();
    vert_handles.resize(3 * num_elems);
    offsets.clear();
    offsets.resize(num_elems + 1);
    for (size_t idx = 0; idx < num_elems; idx++)
    {
        int fid = handle2elem(elem_handles[idx]);
        assert(fid < my_grid->NF);
        assert(is_face_on_surface(fid));

        face_t *f = my_grid->faces + fid;
        int *fnds = f->nodes;
        if (swap[fid])
        {
            vert_handles[3 * idx + 0] = node2handle(fnds[0]);
            vert_handles[3 * idx + 2] = node2handle(fnds[1]);
            vert_handles[3 * idx + 1] = node2handle(fnds[2]);
        }
        else
        {
            vert_handles[3 * idx + 0] = node2handle(fnds[0]);
            vert_handles[3 * idx + 1] = node2handle(fnds[1]);
            vert_handles[3 * idx + 2] = node2handle(fnds[2]);
        }
        offsets[idx] = 3 * idx;
    }
    offsets[num_elems] = 3 * num_elems;
    return;
}

void Mesh_surf3d::elements_get_topologies
(
 const ElementHandle *element_handle_array,
 Mesquite::EntityTopology *element_topologies,
 size_t num_elements, Mesquite::MsqError &err)
{
    (void) element_handle_array;
    for (size_t idx = 0; idx < num_elements; idx++)
        element_topologies[idx] = Mesquite::TRIANGLE;
    return;
}






